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1. Introduction 

Monte Carlo methods for lattice field theories with massless fermions in three or more di- 
mensions continue to pose a variety of challenges. For example, the popular Hybrid Monte Carlo 
(HMC) method encounters many difficulties. Sometimes the method encounters sign problems, 
while in other cases small eigenvalues of the fermion matrix lead to severe singularities. For this 
reason most practical calculations have always relied on extrapolations to the massless limit. Stud- 
ies in Quantum Chromodynamics have shown that reliable extrapolations require calculations at 
many small fermion masses |J^]. As far as we know Monte Carlo calculations on large lattices with 
exactly massless fermions have never been performed so far with traditional Monte Carlo methods 
including the HMC. 

Recently, an alternative Monte Carlo method called the fermion bag approach was proposed to 
solve some four-fermion lattice field theories with exactly massless fermions It was shown that 
the method is extremely efficient at strong couplings where traditional Monte Carlo methods fail. 
Here we argue that the method is also quite general and can be applied to a variety of problems and 
in addition, contains an interesting strong-weak coupling duality which makes the method efficient 
even at weak couplings. The efficiency is due to the fact that the required effort to perform a single 
local update scales like the square of the number of fermion degrees of freedom inside a bag instead 
of the space-time volume. Interestingly, the bag size is a small fraction of the thermodynamic 
volume at both strong and weak couplings. In the massless lattice Thirring model which we study 
here as a first application, the fermion bags typically contain only an eighth of the total degrees of 
freedom even at the quantum critical point. Due to this feature, for the first time we are able solve 
a three dimensional lattice field theory containing exactly massless fermions close to an interesting 
quantum critical point on lattices as large as 40 3 with modest computing resources. Through a 
careful finite size scaling analysis we are able to extract the critical exponents accurately at the 
quantum critical point in this model. 



2. The Fermion Bag Idea 

The idea behind the fermion bag is to identify fermion degrees of freedom that cause sign 
problems and collect them in a bag and sum only over them. This is in contrast to traditional 
approaches where all fermion degrees of freedom in the entire thermodynamic volume are summed 
over in order to solve the sign problem. In some four-fermion models, the degrees of freedom 
inside a fermion bag often involve only a small fraction of the total number of degrees of freedom 
in the entire thermodynamic volume and can be summed up quickly yielding a positive weight. 
This makes the fermion bag approach very powerful. It is an extension of the meron cluster idea 
proposed some time ago [Q]. 

In order to illustrate the idea consider models formulated with n tastes of massless staggered 
fermions. These models contain 2n Grassmann variables per site which will be denoted as Yi( x ) 
and Yi( x ) where i = 1,2, ..n represent taste indices and x denotes the space time lattice point. Let D 
be the V x V free staggered Dirac matrix whose matrix elements are denoted as D^. The properties 
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strong coupling fermion bags Weak coupling fermion bag 



Figure 1: An illustration of a "fermion-bag" configuration at strong couplings (left) and weak couplings 
(right). The interaction sites are represented through sold bonds and the fermion bags are represented through 
the shaded region. At strong couplings the fermion bag are made up of free sites and breaks up into many 
disconnected pieces as seen in the left figure, while at weak couplings the bag contains interaction sites. 



of D are such that any fc-point correlation function of chiral condensates 

Ci(xi,...,x k ) = / [dydyjexp \Y xj/^x) Yi(y))Wi(^i)Yi(^) ■■■ Yi(xk)Vi(xk) (2-1) 

involving the taste i is always positive. Using Wick contractions it is easy to prove that 

Ci(x u ..,x k )=T>et(D) Det(G[{x}]) (2.2) 

where G[{x}] is the k x k matrix of propagators between the k sites Xi,i = l,..,k whose matrix 
elements are G XjrXj = D~^.. It is also possible to argue that 

C i {x 1 ,..,x k )=Det(W[{x}]) (2.3) 

where the matrix W [{x}] is the same as the matrix D except that that the sites {x} = {x,-,i = 1,2, ..k} 
are dropped from the matrix. Thus, W is a (V — k) x (V — k) matrix. The identity 

Det(D) Det(G[{4]) = Det(W [{*}]) (2.4) 

is the basis behind the weak coupling-strong coupling duality we mentioned in the previous section. 

Now consider a generic four fermion lattice field theory action involving n massless staggered 
fermions given by 

S = - £ Wi(x) DxyYi(y) - EE Uij,{xy}Yi(.x)Vi(x)¥j(y)V'j(y) (2-5) 

x,y,i (xy) i.j 

where (xy) refers to some well defined set of neighboring lattice sites. Further we will assume that 
all the couplings Ujj ^ are non-negative real constants. The partition function of the model 

Z= f [dxj/dy] e~ s (2.6) 
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can be expanded in powers of the couplings Uij/^ and each term in the expansion consists of 
a product of correlation functions of the type Q(xi, ...,->Qt) introduced above. Mathematically the 
expansion looks like 

Z= £{[t/]r flQ(x u ...,x kt ) (2.7) 

[{41 '=1 

where refers to some generic power of the coupling and k(,i = 1,2, ..n refers to the number 

interaction vertices for each taste. Note that on a finite lattice the expansion is a very high order 
polynomial and so no convergence issues arise. 

An intuitive physical picture emerges from the above expansion of the partition function. Since 
the kj interaction sites contain both an d Wv fermions are already paired on these sites and do 
not cause sign problems 1 . On the other hand, unpaired fermions that hop freely on the remaining 
sites can indeed cause sign problems. The free sites are collectively referred to as a fermion bag. 
The summation of all fermion world lines inside the bag should be a determinant of a (V — k{) x 
(V —k{) matrix. Indeed Cj(xi,X2, .-,Xjt) = Det(W[{x}]). The determinant can be evaluated easily 
if (V — k,) is small, which naturally occurs at strong couplings. Hence we refer to these bags as 
strong coupling fermion bags. It was shown in [||] that at strong couplings a fermion bag splits 
into many small disconnected pieces making things even simpler. The left figure of Fig. [j] gives 
an illustration of the disconnected pieces of a fermion bag at strong couplings. The solid bonds 
represent interaction sites while the sites in the shaded regions are free sites and form the fermion 
bag. The arrows inside the shaded regions is an illustration of free fermion world lines inside the 
bag. 

Clearly at weak couplings the number of free sites grows enormously. Thus, the definition 
of a fermion bag as a set of free sites, which was natural at strong couplings loses its charm at 
weak couplings. However, thanks to a duality we can construct the fermion bag differently. At 
weak couplings we can view the interactions as the unpaired fermionic degrees of freedom that 
cause fluctuations over the paired fermionic free vacuum as they hop from one interaction site to 
another. The effort to compute these fermionic fluctuations should only grow as the determinant 
of a ki x kj matrix. The duality relation of Eq.( [2.4[ ) shows that Q(x\ ,X2> ■■,Xk) = Det(W[{x}]) = 
Det(D) Det(G[{x}]) where the fluctuation matrix G[{x}] is indeed a x matrix. The right 
figure of Fig|l] gives an illustration of a weak coupling fermion bag. The solid bonds again represent 
interaction sites which form the fermion bag. Fermions hop from one interaction site to another 
through the free fermion propagator. One set of hopping is shown by arrows in the figure. 

We must acknowledge that the weak coupling fermion bag idea is exactly equivalent to the 
idea of summing over all Feynman diagrams and was introduced earlier in the framework of dia- 
grammatic determinantal Monte Carlo methods. For a review please see Ref. [Qj. Indeed the sum 
over all Wick contractions at the given order in perturbation theory is simply Det(G[{x}]). On 
the other hand, we think the fermion bag approach as a more natural interpretation at least in the 
context of lattice field theories since it uncovers the powerful concept of duality discussed above 
and shows that the approach is efficient both at weak and strong couplings. 

The fermion bag approach is rather general and is equally applicable to non-relativistic fermions 
and models with Wilson fermions. However, sometimes it is necessary to introduce unconventional 

1 Any remaining sign problem would also be present in a bosonic field theory. 
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interactions like eight fermion interactions. In some models, the weight of a fermion bag is no 
longer a determinant but involves new mathematical structures like fermionants These are very 
similar to determinants, and may be exponentially difficult to compute [§]. For these models, the 
fermion bag approach is not practical. 



3. Massless Thirring Model 

As a first application of the fermion bag approach we have studied the three dimensional 
massless lattice Thirring model. It is a lattice field theory containing massless fermions and an 
interesting quantum critical point. Variants of this model have been used recently to describe the 
critical points associated with Graphene The model is constructed with two Grassmann valued 
fields y( x ) an d W( x ) on eacn s i te * of a cubic lattice. The lattice action is given by 

s = -Z TO D *y vGO ~ U L TO TO W(y)v(y)- O-i) 

x,y (xy) 

Here (xy) refers to the set of nearest neighbor sites across a bond. We use anti-periodic bound- 
ary conditions in all the three directions. The four fermion coupling U generates the current- 
current coupling of the Thirring model in the continuum. The lattice model is invariant under a 
Uf(l) x U x (l) symmetry, where t//(l) is the fermion number symmetry and U x (l) is the chiral 
symmetry. When the coupling U is small the model contains Nf = 2 flavors of massless four- 
component Dirac fermions at long distances due to fermion doubling. At large U the chiral sym- 
metry breaks spontaneously and generates a single massless Goldstone boson while the fermions 
become massive. There is a quantum critical point U c which separates the phase with massless 
fermions from the phase with massless bosons. 

The model has been studied earlier using mean field techniques [|h, and conventional Monte 



Carlo methods [Q, |10|, |llj, |12(]. In particular properties of the quantum critical point have been 
computed. However, all previous work was done in the presence of a fermion mass and on lattice 
sizes which are not very big compared to the correlation lengths introduced due to the presence of a 
fermion mass. Thus, the analysis may contain uncontrolled systematic errors. A fresh Monte Carlo 
method which works directly in the massless limit should be useful. This is what we accomplish 
here using the fermion bag approach. 

The fermion bag approach for the model was first developed in [||]. However, in the previous 
study only strong coupling bags were used since the concept of duality was not appreciated. Here 
we repeat the study with weak coupling bags and are able to push the study to larger lattice sizes. 



In the fermion bag approach the partition function of the model in Eq.(3. 1 ) can be rewritten as 



Z = Y,U Nb Det(W[n}) (3.2) 

M 

where [n] refers to the configuration of interaction bonds nu^ = 0, 1 and Nb refers to the total 
number of bonds. For a given configuration [n], the free sites form the strong coupling fermion bag 
(see the left figure in Fig. [IJ). However this representation of the partition function is not useful 
when Nb is a small fraction of the volume. For example, at the quantum critical point Nb is about 
an eighth of the lattice volume. On the other hand, one can use the duality relation 

Det(W[n]) = Det(D)[Det 2 (G[«])] (3.3) 
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to simplify the computation when Nb is small. Here D is the free staggered fermion matrix on the 
whole lattice and G[n) is the Nb x Nb free fermion propagator matrix between even lattice sites to 
the odd lattice sites of the bonds. This representation is equivalent to a redefinition of the fermion 
bag as the set of interaction sites (see the right figure of Fig. |]). 

4. Observables 

In order to uncover the properties of the quantum critical point we have measured four observ- 
ables: 

1. The simplest is the average number of bonds (Nb)- Fluctuations in this observable allow us 
to monitor equilibration and autocorrelation times easily. 

2. Since we work with exactly massless fermions the chiral condensate will always be zero. 
However, the chiral condensate susceptibility is nonzero. It is defined as 

X = ^Yi{WxWxWy%)- (4- 1 ) 

x,y 

We expect % to scale as L 2 ^ 11 at the critical point. A related quantity is the bosonic two-point 
function: 

CB(\x-y\) = (WxYxWyYy) (4-2) 
We will define the ratio Cg(L/2 — 1)/Cg(l) as a useful way to track autocorrelations. 

3. Another useful observable that measures the onset of chiral symmetry breaking is the chiral 
winding susceptibility (q\). If we define the conserved chiral charge passing through the 
surface S perpendicular to the direction a as 

(lx)a = E £x n*,a (D^)x,x+a + £ 2e x , (4.3) 

xeS xeS 

where e x = (\) Xl+X2+Xi is the parity of the site x, then (qi) = (|La(^)ce)> is expected to be 
independent of L at the critical point. This allows us to determine the critical point accurately. 

4. Since one of the novel features of the current quantum critical point is the presence of mass- 
less fermions at the critical point, the fermion two-point correlation function is expected to 
show non-trivial scaling. Hence we define 

1 3 

C F (d) = - £ (WxVx+d&) (4.4) 
s a=l 

where x belongs to a site with e (x) = 1 and a is a unit vector along each of the three direc- 
tions. We compute the ratio Rf = Cp (L/2 — 1 ) jCp (1) which is expected to scale as L~( 2+ %) 
at the critical point. 
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Figure 2: An illustration of a zero weight configuration due to the presence of a fermion bag with an odd 
number of free sites. The staggered fermion Dirac operator defined on an odd number of lattice sites has an 
exact zero mode 



5. The Algorithm 

We now discuss the details of the Monte Carlo algorithm we have used to solve the model. An 
important technical problem is that some observables like % get contributions from configurations 
that do not contribute to the partition function. For example, the configuration illustrated in Fig. || 
contributes to %, but has zero weight in the partition function. This is because, in the massless 
limit any strong coupling fermion bag with an odd number of free sites has zero weight. However, 
introducing a y(x) y(x) source in each of the two odd bags makes the number of free sites even and 
thereby changing the configuration weight to be nonzero. Hence instead of computing % (which 
can be done in principle [||]) we redefine our partition function to be 

Z 2 = L / d WV (VxVx W y ¥y) e~ S (5.1) 

x,y J 

which always contains two sources. We then redefine % = I//2 where fi is the fraction of the 
configurations which contain sources one lattice unit apart. The redefined % nas all the same finite 
size scaling properties as the original definition of % close to the quantum critical point. Below we 
will describe the Monte Carlo algorithm for simulating the partition function Zj. The algorithm 
for simulating the real partition function will be the subset of updates in which the sources can be 
ignored. 

The algorithm consists of four updates: (1) Bond creation/destruction, (2) Bond translation, 
(3) Source translation, (4) Worm update. Each update begins with a configuration [n] where 2d 
bond variables n xa = 0, 1 at each site x are completely defined. Note that a = ±1, ±2, • • • ± d 
denote the directions (the negative signs indicate negative directions). If a bond exists between x 
to x + a then n xa = 1. In our notation n x a = n x+ & _ a . The configuration [n] also defines source 
variables m x = 0, 1 such that m x = 1 on only two sites in the entire lattice. On these sites n Xi(X = 0. 
The parity of the site x is given by e x = (l) Jt i+- c 2+-+^. 
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5.1 Bond Creation/Destruction 

Bond creation and destruction is accomplished through a simple Metropolis algorithm. The 
exact update is as follows: 

1. With probability 1/2 we propose to destroy a bond. If Nb = the update stops. Otherwise we 
pick one of bonds randomly with probability l/Nb- We destroy that bond with probability 

'W-'*- <5 - 2) 

where W = Det 2 |g([«]) j and W' = Det 2 |g([«']) j where [n'\ is the new configuration ob- 
tained by destroying the bond. Nf denotes the number of the locations where we would have 
been allowed to create a bond starting from the configuration [n'\. 

2. With probability 1/2 we decide to create a bond. Let Nf to be the number of locations where 
we could create a bond. If Nf = the update stops. Otherwise we pick one of Nf locations 
with probability 1 /Nf. We then create the bond at that location with probability 

NfU W 

P(N h ^N h + l) = -L-- (5.3) 

where W = Det {G(H)} andW' = Det 2 {G(M)} where W] is the new configuration ob- 
tained by creating the new bond. 

If we begin with a configuration with a small number of bonds, since Nf 3> Nb the acceptance 
probability P(Nb — > Nb + 1) is close to one and the algorithm creates bonds efficiently. Once NfU 
is of the order of Nb, the system begins to thermalize and the average number bonds (Nb) fluctuates 
only a little. 

5.2 Bond Translation 

In order to reduce autocorrelation times it is useful to have an update which can move bonds 
from one location to another. This update is again a Metropolis update and is meaningful only 
when Nb / 0. We pick an existing bond at random, destroy it and create a bond at another allowed 
location picked at random with probability 

WW] 

P = (5-4) 
W[n] 

where W'W] and W[n] are the weights of the new configuration and old configurations. 

5.3 Source Translation 

We need to update the locations of the two sources in the configurations that contribute to Z^. 
We again accomplish this using a Metropolis algorithm. We pick one of the sources at random, 
destroy it and create it at an allowed site with the same parity as the site where the source was 
destroyed, with probability 

P = ^l (5.5) 
W[n] 

where W'[n'] and W[n] are the weights of the new and old configurations. 
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5.4 Worm Updates 

Configurations with two sources come in two varieties. One set of configurations contain both 
sources in a single strong coupling fermion bag, while the two sources appear in two different 
strong coupling fermion bags in the other set. The above three updates cannot change between 
these two set of configurations easily and this can lead to a long autocorrelation time. The reason is 
that it is impossible to remove the source from a bag and move it to another bag in one step. Only a 
series of correlated moves can accomplish this task. The problem is most severe in the broken phase 
and perhaps close to the quantum critical point. Fortunately, worm updates of the type discussed 



in Ref. fll3Q completely eliminate this problem. Further, they are extremely fast since they do not 



require the computation of any determinants. We perform two types of worm updates : 

1. Bond Update 

(a) We pick a site x at random and define it as the first site. 

(b) If n Xi(X = 0, the update ends. Otherwise, we break the bond by creating a source at x 
and x + a (i.e., set n xa = 0, m x = 1 and m x+ a) and move to the site y = x + a. 

(c) We then pick a direction /I such that x' = y + fl contains a bond {n x < a i = 1) or x' is the 
first site (x' = x). We pick this direction at random from the available choices. 

(d) If n . ~, = 1, then we break that bond and create a bond that connects the site y and x 
and move the source from y to x + a' (i.e., we set n y ^ = 1, m y = and m x , +( ^, = 1). 
We then redefine y = x' + a' and go back to the previous step. 

(e) If x 1 = x, then we remove the two sources at y and x and create a bond connecting the 
two sites (i.e. set m y = 0, m x = and n y41 = 1). Then the update ends. Note that the 
update begins and ends at the same site and hence is a loop update 

2. Source Update 

(a) We pick one of the two sources at site x. 

(b) We look for all directions /J. that contains a bond. If no such site exists the update ends. 
Otherwise a direction is picked at random from the available choices. 

(c) If the site x + {X contains a bond in the a direction we break this bond and create a new 
bond that connects x and x + fx and move the source at site x to the site x + /t + a (i.e., 
set m x = 0, m x +h+6t = 1 an d n x ^ = 1). The update then ends. 

When the source update is repeated many times it can move the source from one bag to another 
bag. The source update is very cheap and can be repeated thousands of times within a second. 

5.5 Determinant Computation 

Except for the worm updates, all other updates require the computation of the determinant of 
G[n] which is an N/, x Ni, matrix where Ni, is the number of bonds. A change in [n] by a single bond 
or source location typically changes one row and one column of this matrix. The most time con- 
suming part of the algorithm is calculating the determinants before and after the change. Changes 
in [n] can lead to singular matrices quite often. Further, sub-matrices of G[n) can also be singular. 
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There are ways to compute the determinant of a matrix quickly if a single row and column are 



changed Q14[]. However, these tricks involve inverses and must be thought through carefully due to 
possible singularities. 

In our work we use LU decomposition with complete pivoting to compute the determinant. 
However, this approach requires 0(Nf } ) computations. Further, a single sweep requires roughly Nb 
bond updates. Hence, naively the computations necessary to accomplish a single sweep will scale 
like 0(N£). However, as we explain below, we can speed up the update by a factor of Nb, thus 
reducing the number of computations for a single sweep to 0(N^). 

The basic idea is that if one decides apriori that a fixed set of / rows and / columns will be 
updated in the matrix G[n], then while performing the LU decomposition we can organize our 
calculations so that we can reuse a large part of the computation. For example, consider the LU 
decomposition of the block matrix G[n] 



G[n] = 1 (5.6) 




where A, B, C and D are respectively Nb — lxNb — l,Nb — lxl,lxNb — l and I x I matrices. If the 
matrix A is not singular, we can compute the determinant of G[n] as 

DetG[n] = DetA • Det(D - CA _1 fl) (5.7) 

Since A is fixed during the update, the determinant of G[n] due to varying the elements in B,C,D 
requires only 0(lN%) computations. If A is singular, then while computing the LU decomposition 
of A we can isolate the zero mode sector and merge it with the part that is being updated. Since the 
zero mode sector of A is always tiny the above approach works well. 

In order to implement the above idea, we divide the full volume into A^biocks such that every 
bond belongs to a unique block. Each block contains roughly / bonds. We then perform the three 
time consuming updates on the bonds associated to a randomly chosen block j. We choose a basis 



such that the matrix G[n] can be written as a block matrix as shown in Eq. (5.6). In particular 
the matrix A contains propagators only between sites of bonds that do not belong to the block j, 
while the matrix D contains only propagators between sites of the bonds in j. The matrices B and C 
contain propagators between the block j and outside. We then compute the LU decomposition with 
full pivoting on the matrix A. This allows us to isolate any singular part of the matrix A if it exists 
and merge it with the matrices B,C and D. Since the singular part also does not change, it can be 
easily taken into account if necessary. The updates described above can be adapted to each block 
by redefining Nb and Nf as numbers obtained within the block /. By choosing Nu ck ~ 0(\/Nb), 
the number of computations for one sweep scales like 0(N^). We believe our method is similar to 
the "fast-updates" algorithm of the usual determinantal algorithm. 

5.6 Performance of the Algorithm 

We have studied fhermalization times and autocorrelation times for two observables, (Nb) 
and Cb{L/2 — 1), in units of a sweep. In our work, a sweep is defined as V /8 updates of all 



types discussed above . In Fig 5.5, we plot the thermalization time in units of a sweep, for (A 7 /, 



"Since worm updates are cheap we perform several thousand worm updates per sweep. 
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Figure 3: Thermalization time for Nb in terms of sweep using Log plot 




Figure 4: Autocorrelation time for Nf, (Left) and Cb (Right) at \x — y\ = L/2 — 1 in terms of sweep. It shows 
the autocorrelation time does not depend on L even at U c 



for different volumes where we initialize configurations with V /24 randomly chosen bonds. The 
figure shows the system at different volumes can be thermalized within roughly the same number 



of sweeps. In Fig p.5| , we plot the autocorrelation time % in units of a sweep. One typically expects 
T U, where z is the dynamical critical exponent of the algorithm. For most local algorithms 
1 ^ z ^ 2. Many efficient cluster algorithms on the other hand are known to have < z ^ 1. 
Surprisingly, our data shows that z ~ even at U = U c . However, we must emphasize that the time 
spent for each sweep grows as (fV) where / < 1 is some fixed fraction at a given value of U. At 
the critical point / ~ 1/8. Despite this bad scaling with the volume, since / is small we have been 
able to obtain results on lattices as large as 40 3 with modest computing resources. 

6. Results 

One of the main goals of our work is to compute the critical exponents at the quantum critical 



point. These have been calculated earlier in Refs. [9L 10, 11, 12] using the HMC method with the 
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Figure 5: Plots of X^I?^, (q 2 x ) and R f L 2+ ^ as a function of U for L = 12, 16,20,24,28,32,36 and 40. 
The solid lines show the combined fit. Based on the fits we find the critical point to be U = 0.2608(2), and 
three critical exponents are v = 0.85(1), tj = 0.65(1) and rj v — 0.37(1). 



traditional approach in which the four-fermion interaction is converted to a fermion bilinear with 
the help of an auxiliary field. However, these calculations were performed in the presence of a 
quark masses and on lattice volumes that were not very big. The presence of two infrared scales, 
namely the fermion mass and the length of the box can be difficult to take into account in the finite 
size scaling relations which can lead to large systematic errors. In our work, since the fermions 
are exactly massless, the analysis is much simpler and cleaner. We focus on %, (q 2 ,) and Rf in the 
vicinity of U c where we expect the following finite size scaling relations to hold: 



x -'i}-* = ^f k \{u-u c )ri 

k=0 

{q 2 x ) = ^K k \{U-U L )L^ 

k=0 

3 r , 

Rfl }+v ¥ = £ Pk (U-U C )L* 



k=0 



(6.1a) 
(6.1b) 
(6.1c) 



For each observable, when the left hand side is plotted as a function of U, curves for different 



values of L must cross at U = U c . A combined fit to all the Eqs. (6.1) gives the critical exponents 
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V = 0.85(1), tj = 0.65(1) and T) w = 0.37(1) and U = 0.2608(2) with a X 2 /d.o.f. = 1.3. The 
complete list of fit parameters are listed in Table [l[ The results from the combined fit are plotted in 
Fig. |5[ For comparison, one of the earlier work finds v = 0.71(4) and v = 0.60(2) JlTjl. 



Table 1: Results for the fit parameters from the combined fit of the data to Eqs. (6.1). 





% 


V 


U c 


Kb 




Kz 


*3 


0.65(1) 


0.37(1) 


0.85(1) 


0.2608(2) 


0.369(3) 


0.63(1) 


0.52(2) 


0.09(1) 


fo 


fx 


h 


h 


Po 


Pi 


P2 


P3 


2.52(3) 


-2.53(5) 


0.71(3) 


0.10(1) 


33.9(2) 


-5.0(1) 


-2.0(2) 


-2.5(5) 



7. Conclusions 

In this work we have shown that the recently proposed fermion bag approach is a powerful 
technique for solving some four-fermion lattice field theories. Due to an interesting duality, the 
approach is efficient both at weak and strong couplings and continues to perform well at intermedi- 
ate couplings. As a first application of the method we studied the critical behavior in the massless 
Thirring model and found an algorithm that practically eliminates critical slowing down when times 
are measured in units of sweep. The time to perform a sweep scales as 0(N^) where Nb is the size 
of a fermion bag which is usually a small fraction of the volume. In the massless thirring model we 
found that Nb ~ V /8 at the quantum critical point. The Hybrid Monte Carlo method has never been 
successfully applied to massless fermions on large volumes. We have accomplished this using the 
fermion bag approach on lattices as large as 40 3 with moderate computing resources. We were able 
to compute the critical exponents at the quantum critical point. 
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